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SUMMARY 


A new implicit algorithm for time step integration of finite element 
structural dynamic equations is presented. Convergence, stability and numerical 
damping properties are discussed. Due to the way nonlinear structural behavior 
is taken into account the algorithm is expected to compare favourably with ex- 
isting ones. Some simple numerical results are presented. A related explicit 
algorithm is also derived and shortly discussed. 


INTRODUCTION 


Lately much attention is being paid to finite element structural dynamic 
problems where nonlinear behavior is taken into account. It is typical for this 
kind of problems that the internal nodal forces developed by the structure and 
resisting external loads and inertia forces are not linear functions of nodal 
displacements but have to be evaluated from the actual stress state of the de- 
formed structure by virtual work integration or similar procedures. The dis- 
placement time history is then obtained integrating step-by-step the nonlinear 
dynamic equations, a very cumbersome procedure in most practical cases. 

Time step algorithms can be classified as implicit, requiring the solu- 
tion of a system of coupled and generally nonlinear equations at each time step 
or as explicit where the unknown problem parameters at the end of each step are 
obtained directly. A clear and concise discussion of requirements and applica- 
bility for both kinds of algorithms can be found in reference [l ] . 

In the present paper an implicit algorithm is first derived and discussed 
using simple numerical tests to show some of its properties. A corresponding 
explicit algorithm [which has not been implemented yet] is then also derived. 

It should be clear, however, that the present paper has the limited scope of 
reporting some early results of a research project presently in progress. More 
extensive numerical tests and comparisons with different algorithms are needed 
for general conclusions. 
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IMPLICIT METHOD 


Within a time step At, i.e. between t D and t ^ = t Q + At (with T = t-t Q ), 
the following system of differential equations, a dot indicating derivation with 
respect to time, has to be solved: 

{R(x)} = [M]{W(x)} + [C]{W(T)} + {F(x)} - {P(x) } = 0 

where: {R(x)} = residual vector; {W(x)} = vector of nodal displacement para- 
meters; [M] = mass matrix; [C] = viscous damping matrix; {F(x)} = vector of in- 
ternal nodal forces; {P(x)} = vector of prescribed external loads. The internal 
structural nodal forces (F(x)} are evaluated as follows (see Fig. 1): 

{F(t)> = {F } + [K] ( {W(x) } - {W }) + {F(x)} 
o o 

where: {F } = internal forces at t = t Q (subscripts 0 and 1 always refer to the 
beginning and the end of the time step); [K] = an approximation of the actual 
secant stiffness matrix between t Q and t^ (generally the tangent stiffness matrix 

evaluated at t = t Q or at the begin- 
ning of some earlier time step will be 
used for [K.]);{F(x)} = vector of cor- 
rective internal forces due to the 
nonlinear behavior of the structure 
within the time step. The distribu- 
tion of these corrective forces {F(t)} 
is assumed to be given by: 

{F(x)} = fb ii (x)J{F 1 } 

where the b-^Cx), coefficients of the 
diagonal interpolation matrix Hb.j^(l)3, 
are time functions (with b-j_j_(0) = 0 
and bj_i(At) =1) to be chosen according 
to the expected distribution of each 
single corrective force Fi(x) as ex- 
plained later. The internal corrective 
forces {F-]} at the end of the time step 
will have to be evaluated by stress-virtual-strain integration or by similar 
procedures (taking into account strain-stress history, large displacements, etc.) 
for all elements where nonlinear behavior is expected. 

For the displacement parameters the following assumptions are used: 

{W(x ) } = N (X)*{W } + IM^ (X ) * ( W q } + N (x)-{W } 

the Ni’s being shape functions in the time dimension. In order to insure the 
necessary continuity of {w} and {w} (but not of {w}) between time steps the fol- 
lowing conditions have to be satisfied: 



142 



N (0 ) = 1 
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We also impose the following conditions: 
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At 
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the parameter p CO < p _< 1 ) being the spectral radius in the limit At -* 03 of the 
operator matrix for the linear case as explained later. Polynomials of 4th degree 
in s (with s = x/At) can be chosen for the shape functions Nj_: 


l\l 2 = C-5C1+p]s‘ t + 4 ( 3 + 2p ] s 3 - 3 ( 3+p ) s 2 + 2s)-~ 

N = ( 5 C p 2 -4p+7 3 s 4 - 8 ( p 2 - 5p+9 ) s 3 + 3 C p 2 - 6p+1 3 ] s 2 ) •— - 
o Z Z p 


For the important special case where p 
of the Ni functions is shown in Fig. 2. 


1 (no numerical damping) the shape 


N,(r) N 2 (t) 



N 3 (r) 

\ 

I 
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Fig. 2 : Shape Functions N, (r) N 2 (r) N 3 (r) 


In order to evaluate { } = { W ( t ^ ) } the weighted integral of the residual 
vector {RCt ) } is set to zero at each time step: 


At 


Jg* { R ( x ) } • dx = 0 


with: G = 


At 

1 +p 


where the choice of a constant weighting function G insures that the integral of 
all external and internal forces vanishes at each time step. The following system 
of equations for {AW} = { } - {W Q } is found: 

C [M3 + ^-[C] ♦ 


1 +P 


1 ' 

T i^ [R,H4W > ■ At[MH “o > - - < p *» 


where : 


At 


fB. J = — - Jtb. . (x)J*dx 
ii At J li 


At 


{P*}= ~ /{P(x)}*dx 
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As the matrix multiplying the unknown vector {AW} has the same structure 
as the stiffness matrix [K] , the algorithm is called implicit. The corrective 
forces {F^} can only be evaluated when {AW} is known; the system of equations is 
therefore nonlinear. Unless the influence of {F^ } is neglected, equilibrium it- 
erations within each time step are necessary, requiring reiterate evaluations of 
{F-]}. A norm for the changes of {F^} (rather than {Aw}) was used in all numeri- 
cal examples as a convergence criterion. 

The displacement velocity vector { W-^ } is found from the assumed shape func- 
tions: 

{ V = 1? {AW} " p{ V 

The convergence and stability properties of the algorithm for proportional- 
ly damped, linear-elastic free vibrations can be discussed applying a modal trans- 
formation (see reference [2]). For a mode Y(t) associated with a frequency to and 
with a damping coefficient E, the following operator relations are obtained 

(ft = w*At ) : 



The operator matrix shows the following properties: 


1. For small At ' s i.e. in the limit ft -* 0 convergence to the true solution is 
insured . 

2. The factor p represents the spectral radius of the two-roots operator matrix 
in the limit ft a- oo. The algorithm is therefore unconditionally stable for 

p j< 1 . 

3. With p = 1 the operator matrix becomes identical to Newmark's with y = 1/2 and 
6 = V 4 ("trapezoidal rule”, see references [3] and [2]]. No numerical damping, 
second order accuracy, minimal period elongation and no amplitude decay are 
obtained . 

4. For p < 1 numerical damping used to "filter out” disturbing high frequencies 
is introduced. In fact, the choice of p represent a convenient and natural 
way of controlling numerical damping. This, however, decreases linearly with 
ft (as in Newmark's method with y > V2) and not quadratically as might be 
desirable (see references [2], [4], [5]). For p < 1 second order accuracy is 
therefore lost, which limits the choice of p to values close to unity. 

5. The so-called "overshoot" effect both for displacements and velocities (see 
reference [6]) is avoided for all values of p. 

All this is of course only valid in the linear case. However, the algorithm 
has two other properties worth mentioning as they appear to be valuable in the 
nonlinear case. 
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The first one concerns the way the b^-jjx) functions describing the dis- 
tribution of the internal corrective forces {F(x)} within the time step or the 
corresponding average values Bn are chosen. 

By setting B.^ = 0 the influence of the corrective force F^ is neglected, 
i.e. linear behavior as described by the stiffness matrix [ K ] is assumed in that 
zone of the structure which affects F^. If all Bj_^ are set to zero, { F^ } is neg- 
lected and no equilibrium iterations are necessary, which may lead to rapidly 
diverging results unless sufficiently small time steps are used. 

If nonlinear behavior is expected and the tangent stiffness matrix evaluat- 
ed at the beginning of the time step is used for [ K ] then Bj_i = 1/3, correspond- 
ing to bj_j_ = (x/At) 2 , should be used, because the time derivatives Fj_[0) of the 
corrective functions at X = 0 are Known to vanish. 

If for [K] the tangent stiffness matrix of an earlier time step is used, 
then Bj_j_ = V2 may be used corresponding to the assumption of a linearly distri- 
buted corrective force T ± , i.e. to b i± = x/At. In the commonly used undamped 
version of the Newmark algorithm (with y = 1/2 and 3 = 1/4) implemented with the 
so-called "out-of-balance-load" procedure as well as in the original iterative 
formulation of the Newmark method (see reference [3]) the Bj_j_’s are implicitly 
always set to V2, a rather poor choice when the actual tangent stiffness matrix 
is used. 

In some parts of the structure, where highly, nonlinear behavior is expect- 
ed, it might be usefull to evaluate the F^'s not only at x = At but also at 
x = At/2 or even at several points between x = 0 and x = At . This allows evalu- 
ation of the Bii coefficients more exactly by simple numerical integration pro- 
cedures . 

The second welcome property of the algorithm is due to the fact that for 
each time step average values {P*} of the external loads {P(t)} are used in- 
stead of considering only instantaneous values of { P ( 1 3 } as in most other well- 
known algorithms. If the time step At is large compared to the typical period of 
the loads and if the load values are known in intervals smaller than At, quite 
relevant improvements can be obtained by using very little-time-consuming nu- 
merical integration procedures to evaluate {P*}. 

NUMERICAL RESULTS 


The single-degree-of-f reedom, numerically and physically undamped (p =1; 

E, = 0) examples presented below are intended to show the beneficial influence of 
choosing the proper Bj_j_ coefficients and of correctly evaluating the average ex- 
ternal loads {P*}. 

In the first two examples the homogeneous equation of motion is given by: 
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M«W + F(W] = 0 


with a constant mass M, a prescribed starting velocity W(o] and a vanishing dis 
placement W(o] at t = 0. In the first example a quadratic relation between F(W) 


and W (nonlinear elasticity] is assumed 
pie the internal force F(W] is produced 



Fig. 3 : First Example 

F(w)-w-Relation 


as shown in Fig. 3. In the second exam- 
by five linear-elastic brittle springs 



Fig 4 : Second Example 

F(w)-w- Relation 


with identical stiffnesses. Four of them are assumed to break at W = 6, 26, 36 
and 46 leading to the F-W-diagram shown in Fig. 4 remindful of crack propagation 
problems . 

Different time histories for different B..’s are shown in Figs. 5 and G. 

li 
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1 2 3 


Fig. 6 : Second Example. Time History for Several B u Coefficients 

The "exact" ones were obtained using very small time steps. For the others large 
time steps were used. The tangent stiffness Ci.e. the derivative of F(W] with 
respect to W) was updated at each time step and equilibrium iterations [except 
for = 0] were repeated until the change of the corrective force F| became 
negligible. "N" is the average number of such equilibrium iterations. "FI" is 
the number of times the internal corrective forces {"F(t)} have to be evaluated 
at each time step for the different B^'s: F! = I\1 = 0 for = 0 ; FI = 1 for 

B ii ~ 1/3 or B ii = 1/2: for FI > 1,B-j_i is variable as FI values of {F(t)} are used 
to evaluate the time integral leading to Bj_j_. 

The advantage, when using the actual tangent stiffness matrix, of choosing 
Bj_i = 1/3 instead of B-^ = 1/2 (which, as explained earlier, would correspond 
to the most used version of Newmark's algorithm) is evident. An even more marked 
improvement, as clearly demonstrated in Fig. 6, is obtained by evaluating Bj_i 
more exactly [FI > 1). This, however, requires time-consuming evaluations of 
iF(T)} for different t’s within each time step. 

The third example (Fig. 7) shows the response of a linear undamped system 
with a natural period T s = 4 to a sinusoidal load with a period Tp = 1 . The time 
step used (At = 0.3) is very large compared with T but reasonable if compared 
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with T s . The parameter ”K" given in Fig. 7 represents the number of discrete 
values of the load function P(t] used for the evaluation of the load average P* 
within each time step (K = 2 if only P C 1 0 3 and P C t ^ 3 are used as in most well- 
known algorithms]. Of course convergence to the exact solution due to the large 
time step used is not obtained; nevertheless the advantage of correctly evaluat- 
ing P* is evident, the computational effort needed being negligible. 



o 


EXPLICIT METHOD 

An explicit algorithm closely related to the implicit one discussed above 
has also been derived. The main difference lies in the way the integral of the 
internal forces { F C T 3 } within the time step is evaluated: 

At 

/{F(x ] }*dt 2L At * (f) 
o 

where {F} is determined assuming a linear distribution of (F[t)} as well as 
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uniform motion within the time step: 


{F} = { F C { W} 3 } with: {W} = {W } + ~{W } 

o 2 o 

Requiring the weighted integral of the residuals to vanish while using the 
same shape functions and the same constant weighting function G as before, 

the following system of equations for {AW} is obtained: 

C[M ] + ^[C]){AW} = At[M]{W Q } - f^({F} - {P*>3 

Assuming the matrices [M] and [C] to be diagonal, i.e. lumping masses and 
viscous dampers in the joints, these equations are trivial. The algorithm is 
then called explicit. 


If, as for the implicit method, a modal transformation is performed, the 
following operator equations are found: 



At * 


1 

1 + T 2 — 
1 +P 


1 , Q 2 

1+p 1 +p 

-Q . 2 


2fflp _ 

1 +p 2 


< 


Y 

o 


> 


At • Y 

o 


As expected the operator matrix shows that the algorithm is only condi- 
tionally stable. For p = 1 the stability condition obtained by setting the 
spectral radius equal to unity is: 

=2 

cr 

or, for the more stringent "bifurcation” condition (see reference [ 7 ] 1 : 


n 


bif 




For damped systems ( E, > 0] these conditions are less stringent than those 
given in references [7] or [8], due to the fact that here a diagonal damping 
maxtrix [C] is assumed. For 0.8 <_ p < 1.0 the changes in f2 cr and i2bif are found 
to be minimal. However, as in the implicit method, with p < 1 numerical damping 
is introduced and second order accuracy is lost. The parameter p has not the 
same meaning as before but still controls numerical damping. 

In fact it is questionable if a choice different from p = 1 would make 
much sense when using the explicit method. Because, due to stability reasons, 
explicit methods require very small time steps, they are mainly used for problems 
where all or almost all frequencies of the finite element model have to be taken 
into account (e.g. for shock or wave propagation problems] so that a "filtering 
out" of frequencies by numerical damping is not necessary and mostly not desir- 
able. 
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The main reason why p was left as a variable in the above derivations is 
due to the fact .that it seems possible to combine the implicit and the explicit 
method following the procedure suggested by Hughes et . al . in reference [8], 
where their damping coefficient y is also used in both cases. 


FURTHER DEVELOPMENTS 

As stated in the introduction this paper is only concerned with early 
results of a research project presently in progress. The following developments 
are planned: 

1. Extensive tests of the implicit algorithm for reinforced concrete structures 
under earthquake loads. Different element formulations are to be compared. 

2. Implementation and tests of iterative methods (e.g. overrelaxation] for the 
solution of the equations arising at each time step, possibly combining iter- 
ative solution steps with equilibrium iterations. As good guesses for start- 
ing iteration are available from earlier time steps, iterative methods might 
well prove to be quite efficient. 

3. Implementation and tests of the explicit algorithm and of the coupled im- 
plicit-explicit-method as suggested by Belytschko and Mullen in reference 

[ 8 ] or by Hughes et . al. in reference [7j. Also in this case the use of iter- 
ative methods for the solution of the implicit equations looks promising. 

4. Development of efficient, state-of-the art computer programs to be applied 
in practical cases. This is a critical point, the applicability of step-by- 
step procedures being severely limited by the great amount of computations 
involved. 
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